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We present a quantum Monte Carlo study of the hydrogen-benzene system where binding is 
very weak. We demonstrate that the binding is well described at both variational Monte Carlo 
(VMC) and diffusion Monte Carlo (DMC) levels by a Jastrow correlated single determinant gem- 
inal wave function with an optimized compact basis set that includes diffuse orbitals. Agreement 
between VMC and fixed-node DMC binding energies is found to be within 0.18 mHa, suggesting 
the calculations are well-converged with respect to the basis. Essentially the same binding is also 
found in independent DMC calculations using a different trial wave function of a more conventional 
Slater-Jastrow form, supporting our conclusion that the binding energy is accurate and includes all 
effects of correlation. We compare with empirical models and previous calculations, and we dis- 
cuss the physical mechanisms of the interaction, the role of diffuse basis functions, and the charge 
redistribution in the bond. 

PACS numbers; 31.15.ae, 02.70.Ss, 68.43.Bc 



I. INTRODUCTION 

A great deal of research has gone into the study of hy- 
drogen storage materials due in large part to the prospect 
of zero emissions transportation. The value of find- 
ing a material that effectively and reversibly stores hy- 
drogen can hardly be overstated, [ij Hydrogen storage 
has a long history starting with the discovery of re- 
versible hydriding in palladium in 1866 by Graham, % 
with much work on systems such as metal hydrides 3j 
and complex hydrides. 0, 0] Recently, materials 

such as carbon nanotubes, M fullerenes, [§| metal-organic- 
frameworks (M0F),[l3, [m, and others have been 
studied. The U.S. Department of Energy has set a goal 
for 2010 to find a material that can meet the requirements 
of storing molecular hydrogen at 6% weight and 45 grams 
per liter, which is nearly half the hydrogen content of wa- 
ter, reversibly in the range of -30-50°C for a thousand 
cycles, [ijli While many solutions to this problem have 
been offered, none have satisfied all these constraints. 

In this paper we study the adsorption of molecular 
hydrogen on a benzene ring as a test system. By it- 
self this system is not expected to be a practical system 
for storage due the expected binding energy being much 
weaker than the target range of 20-40 kJ/mol H2 7-15 
mHa/H2) needed for reversibility in the desired temper- 
ature range, llj However, benzene is similar to the 5- or 
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6-member rings that are characteristic building blocks of 
all the carbon systems mentioned above. As such, an 
accurate description of this structure is highly relevant 
to ongoing research where hydrogen is adsorbed on or 
around carbon rings. Besides this, the hydrogen-benzene 
system is a good test case for theoretical predictions be- 
cause of the stringent requirements to reliably determine 
binding energies at the desired accuracy level. Further, 
a careful study of this system is an im por tant test of the 
transferability of empirical potentials|l4l. fisl l that have 
been constructed primarily from experimental data on 
graphitic systems. In this paper we consider only the case 
where the hydrogen dimer is oriented along the Cssym- 
metry axis of the benzene molecule. Other papers [l6l[l7| 
have found this to be the favored configuration and ori- 
entational differences are not taken into account in this 
work, since our main purpose is to present benchmark 
calculations for the most stable geometry. 

There have been many previous studies of the bind- 
ing of H2 on benzene and related systems using vari- 
ous methods including density functional theory (DFT), 
M0ller-Plesset second order perturbation theory (MP2), 
coupled cluster (CC) with single and double excitations 
(CCSD), and variations of these. [H [H, 113 The values 
obtained so far for the binding energy, falling in the range 
of 0.4-1.9 mHa, [ll|, El, [I3|, are very small and require 
a high level of accuracy of all the methods. The DFT 
calculations have great advantage as they are fast, scale 
well with system size, and can be readily converged with 
respect to the basis. However, the accuracy of their re- 
sults is limited by the approximation on the exchange 
and correlation functionals, and there is no known way 
to systematically improve it. The many-body CC meth- 
ods are the most accurate, although their applications 
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are limited to small systems and not-so-large basis set 
due to poor scaling with the number of electrons and 
the size of the basis. Perturbation methods such as MP2 
are valuable theories, with a system size scaling better 
than any CC method, but with intermediate accuracy. 
In the H2-benzene system, perhaps the most accurate re- 
sults to date have been derived from MP2 and CCSD(T) 
calculations by Hiibner et o/.[l3] They found that the 
binding increases with increasing basis size in MP2 cal- 
culations, whereas it decreases as the level of the theory is 
improved to CCSD(T). Based on the best MP2 binding 
energy (1.87 mHa) and the best CCSD(T) value (1.16 
mHa) with affordable bases, they estimated the actual 
binding energy to be ~ 1.5 mHa. It should be noted 
that these numbers have already included basis set su- 
perposition error corrections as high as 0.36 mHa and 
so represent a substantial fraction of the binding energy. 
Such methods necessitate carefully extrapolating the re- 
sults with respect to the basis set and the level of theory. 
However, such extrapolations represent a difficulty due to 
the high computational cost of large bases, particularly 
in the CCSD(T) framework. 

In the present work we study the hydrogen-benzene 
problem using quantum Monte Carlo (QMC) methods, 
that offer several advantages: many-body correlation ef- 
fects can be explicitly included in the wave function, scal- 
ing with the system size is favorable like DFT, and cal- 
culations are variational and usually less dependent on 
the basis set. For a review and references to earlier work, 
see Ref. [l^. By means of QMC, a trial correlated wave 
function can be optimized in the variational Monte Carlo 
(VMC) framework, [H HO] and its energy can be further 
minimized by the diffusion Monte Carlo (DMC) algo- 
rithm, which stochastically projects the optimized VMC 
(trial) wave function to the ground state. The only fun- 
damental limitation is the well known "sign problem" for 
fermion systems, that does not allow a numerically sta- 
ble calculation. Therefore, in this case, the so called fixed 
node (FN) approximation is adopted by constraining the 
diffusion within the nodal pockets of the initial varia- 
tional wave function. [11] Thus, the FN DMC method is 
unbiased only if the nodes of the trial wave function co- 
incide with those of the true ground state. We addressed 
the issue of the FN bias in two ways. First, we used ad- 
vanced QMC optimization methods [Tol [20l . [2ll . [2^ and 
physical principles to find a Jastrow correlated antisym- 
metrized geminal product (JAGP)|2l, [2^ which gives a 
VMC binding energy with an accuracy comparable to 
the post Hartree-Fock (HF) methods. In addition, we 
computed the binding energy at the DMC level using 
the JAGP and a simpler Slater- Jastrow (SJ) form with a 
PBE-DFT optimized basis set as the trial wave function. 
The agreement found between them supports clearly the 
idea that our results are independent of the basis set and 
variational form, and it is a check for the accuracy of our 
DMC calculations against the FN approximation, since 
the nodes of the two wave functions are a priori different. 
This is encouraging for another reason: although the SJ 



trial function is not as accurate as the JAGP at the VMC 
level, it is more easily extended to larger systems which 
are important for future work. The necessary condition 
for that agreement is using a basis sufficiently extended 
in the tails. This is not surprising for a system driven by 
Van der Waals (VdW) interactions which lead to weak 
binding and large equilibrium distance, as pointed out in 
Ref. 25, but it is a crucial point since the tails are not 
very important in the total energy. Our work shows that 
DMC can capture the correct binding as long as the ba- 
sis is extended enough to allow accurate sampling of the 
outer regions of the molecules. This is brought out by 
a detailed study of the electron density changes due to 
binding. 

The paper is organized as follows. In Sec. |TT] we de- 
scribe the QMC methods employed as well as the SJ and 
JAGP wave functions that serve as the variational guess 
in our QMC calculations. In Sec. IIIII we discuss our re- 
sults on the binding energy of the hydrogen-benzene sys- 
tem, where the hydrogen is oriented perpendicular to and 
centered over the benzene at various molecular spacings. 
In Sec. IIVI we compare our findings to previous works. 
In Sec.|V]we discuss the physics of the hydrogen-benzene 
bond in terms of its electron density, by comparing the 
QMC and DFT-PBE results. Finally, we draw our con- 
clusions in Sec. IVII 

II. COMPUTATIONAL DETAILS 

The unique feature of QMC methods is that they in- 
volve directly the many-body wave function and address 
the necessary high-dimensional integrals by stochasti- 
cally sampling the configurational space. [H, [2^ In so 
doing, the wave function may take on a form more gen- 
eral than a bare linear combination of Slater determi- 
nants while still maintaining computational efficiency. In 
the following Subsections we describe the variational trial 
wave function (Subsec. Ill Ap . and provide a brief intro- 
duction on the methodology (Subsec. Ill B[) . 

A. Wave functions 

A wave function that describes a system of A'' identical 
fcrmions must be antisymmetric under particle exchange. 
To simplify the description of such a wave function, it is 
often useful to factor the wave function into a positive 
symmetric part, called the Jastrow factor, and an anti- 
symmetric part so that a wave function can be expressed 
as 

*(xi, . . . ,XAr) = J(xi, . . . ,XAr) *As(xi, . . ■,Xn)- (1) 

where x^ = {ri,ai} is a space-spin coordi- 
nate, J(xi, . . . , x^v) is the Jastrow factor, and 
^^^(xi, . . . , xjv) is the antisymmetric part. The Jastrow 
can be further factored into one-body, two-body, three- 
body, and higher-body terms {J = JiJ2«/3''') which 
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correspond to effective electron-ion, electron-electron, 
electron-electron-ion, etc. interactions. 

One of the choices of trial function is to approximate 
the antisymmetric wave function as a single Slater deter- 
minant of spin orbitals. If there are no spin orbit inter- 



J 



actions, the energy depends only upon the spatial part of 
the wave function which can be written as a product of 
spin up and spin down determinants. In the unpolarized 
case, the spatial form is given by 



5'As(xi, . . . ,XAr) 



I 



V"/2+i(ri) ... ipN{r{) 
'^«/2+i(ri/J ... ^Niri^J 



(2) 



where each (pi{r) is a single-body space orbital and N 
is the total number of electrons. In our calculations the 
single-body orbitals u^j (r) are derived using the Perdew- 
Burke-Ernzerhof[23, (PBE) density functional with 
the VTZ Gaussian basis [2^] modified to include diffuse 
functions from the aug-cc-pVTZ basis, (soj The normal- 
ization factor in Eq. [2] has been dropped since in the 
QMC approach the wave function does not need to be 
normalized. 

Since the Slater determinant description alone is un- 
able to treat correlation effects, the Jastrow factor is used 
to include electron correlation in a computationally effi- 
cient way. Moreover, it provides a way to enforce cusp 
(or coalescence) conditions which are an exact property 
of the many-body wave function, and demand that the 
local energy El = '^^^H^> remains finite as electrons 
approach ion centers and each other. The Jastrow factor 
we applied to the Slater determinant is a Wagner-Mitas 
form _31] modified so that the electron-ion and electron- 
electron cusp conditions are fulfilled. The one- and two- 
body Jastrow terms are given by 



Ji(R) ==[]exp 



'^{bakna + Cak)Vak{na) 



and 



J2(R) = nexp 

i<j 



E 

. k 



[hkTij + Ck)vk{rij) 



(3) 



(4) 



where R — {ri, . . . , rjv} specifies the N electron space co- 
ordinates, i and a index electrons and nuclei respectively, 
ria and are electron-ion and electron-electron dis- 
tances, and k indexes the expansion terms. In our work 
we used three terms and, when needed, one cusp term. 
In the above Equations, Vk{r) = (1 — z(r/rc„t))/(l -f 
(3kz{r /rcut)), with z(x) = x^(6 — 8a; + ix^) and parame- 
ters h, c, (3 optimizable (with the exception of those that 
are cusp dependent). The function z{x) has the prop- 
erties z(0) = z'(0) = z'(l) = and z{\) = 1, so that 
the Jastrow has a well defined cutoff at rcut = 10 Bohr. 
Cusps between same spin electrons are not accounted for. 



This is justified because of the Pauli exclusion principle, 
which keeps them apart. Also the three-body terms have 
been neglected. It should be emphasized that the single- 
body Slater orbitals obtained from PBE-DFT are not 
further optimized since we would like to check the ac- 
curacy of the PBE-DFT nodes with respect to a more 
correlated and fully optimized wave function, such as the 
JAGP form described below. However, optimizing the 
above Jastrow is convenient as it improves the VMC en- 
ergy and variance and shortens the DMC projection time, 
without changing the nodes. 

The other trial function used in this work is the JAGP, 
where the antisymmetric part is a single determinant of 
two-body orbitals (geminals). This approach has been 
successfully applied in several contexts where electron 
correlations play a significant role. For example, the 
JAGP form is related to the pairing in the BCS wave 
function for superconductivity, [SZ, 33] the resonating va- 
lence bond (RVB) proposed by Pauhng in 1939,^34] and 
can be used to describe strongly-correlated electrons in 
transition metals. Recent applications in quantum chem- 
istry include benzene. (23j the benzene dimer interacting 
via weak van der Waals forces, [20| and iron dimer. [35| 

Since the ground state of the hydrogen-benzene system 
is an unpolarized spin singlet [N^ = = ^/2) the 
spatial part of the AGP wave function can be written as 
a determinant of pairing functions |,36i| without including 
unpaired orbitals, namely 



'i'As(X) = 



(rl,ri) 



(r3v/.,ri) 



^ N/2-1 ^ N/2) 



(5) 



where X ~ {xi, . . . ,X7v} specifies the N electron space- 
spin coordinates and the paring function (f){v\ , rj) can be 
expanded in single-body atomic orbitals so that 

</.(rJ,r))= 5^ Ai7(^,KrI)^.™(rj) (6) 

Imab 

where I and m index the orbitals centered on ions a and 
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b respectively. Here also, as in the SJ wave function, 
Gaussian type orbitals are used. 

The Jastrow factor in the JAGP wave function is some- 
what different from the one applied to the Slater deter- 
minant. The cusp conditions are fulfilled through the 
one-body Ji, and two-body J2 Jastrow terms, written as 

Ji(R) = Hexp \-{2ZafMi'2Zay^'r,a)\ , (7) 



and 

J2iR)^lle^v[uin,))], (8) 

i<j 

where R is an all-electron configuration, i and j are elec- 
tron indices, and a is a nuclear index. The ion centers 
have effective charge Za and the function u(x) satisfies 
the electron-ion and electron-electron cusp conditions be- 
tween unlike-spin particles with u(0) = and u'(0) = ^. 
Here, u{r) = -j (l — e~'^/^), where F is an optimizable 
parameter. In Eq. [71 the argument of u is multiplied by 
{2ZaY^^ in order to satisfy the random phase approxima- 



tion behavior at large r,v,.|37| 



A distinguishing feature of the JAGP with respect to 
the simple SJ wave function is the presence of electron- 
electron-ion and electron-ion-electron-ion terms, conven- 
tionally referred to as three- and four-body Jastrow fac- 
tors. In the JAGP wave function, they are written as the 
exponential of a pairing function like the one in Eq. [6l 
namely 



J34(R) = J|exp 



ablm 



(9) 



Here (7°^ are optimizable parameters and I (rn) is an in- 
dex for single-particle Gaussian orbitals Xai centered on 
nucleus a (b). The three- and four-body Jastrow terms 
provide for electron-correlations substantially beyond the 
largely cusp related one- and two-body terms and are able 
to describe subtle effects like van der Waals forces, [ssj 
However, Eq.[n]does not include the three-body cusp con- 
ditions recently derived by Fournias et al.,^^ which can 
improve the quality of the nodes of the JAGP wave func- 
tion described here. The effect of the three-body cusp 
conditions in the energy optimization and nodal struc- 
ture is presently under investigation. 



work). We proceed to optimize its energy and variance 
at the VMC level usin g m inimization methods suitable 
for the particular form jl9l [20I |40| . [4l| The resulting an- 
alytic wave function is projected to the FN ground state 
using DMC methods [42, |43| recently developed to yield 
a stable simulation and an upper bound of the ground 
state energy even for non-local pseudopotentials. 

As we mentioned above, we use the full electron- 
nucleus Hamiltonian except for the carbon core which 
is replaced by a pseudopotential. This leads to a bet- 
ter statistics due to a narrower energy scale, a reduction 
in the number of optimization parameters, a more stable 
optimization of our JAGP wave function, [20| and a larger 
DMC time step needed for convergence, which results in 
a cheaper computational cost of the simulation. On the 
other hand, its drawback is that part of the fully local 
Coulomb potential is replaced by a non-local pseudopo- 
tential Kion-iocai that is angular momentum dependent. 
Within the VMC framework the corresponding angular 
integration of the non-local potential remains possible 
since the wave function is known analytically. However, 
problems arise in the FN DMC because the FN ground 
state is given only by a stochastic sampling. A partial so- 
lution is the localization approximation (LA), where the 
trial (or guiding) wave function \E'g is used to approx- 
imate the projected ground state so that the non-local 
pseudopotential terms can be evaluated, [isj 

This changes the Hamiltonian so that the projected en- 
ergy is no longer a variational upper bound of the origi- 
nal non-local FN Hamiltonian. The FN Green's function 
positive definite character becomes impaired by the very 
attractive parts of the non-local pseudopotential which 
has been made local using the guiding wave function. 
Indeed, on the nodes of the guidance the localized poten- 
tial can diverge negatively leading to possible numerical 
instabilities where the walker population can grow with- 
out bound. Our method of implementing non-local pseu- 
dopotentials in the FN DMC methods sidesteps these 
problems entirely. 

Our FN DMC calculations are done with either con- 
tinuous or lattice regularized (LRDMC) moves both of 
which utilize a common means of addressing the inher- 
ent problems of the LA. In contrast to the LA, we use 
a breakup [4^, [43j of the non-local potential that local- 
izes the positive matrix elements into the branching term 
while treating the negative matrix elements as a non-local 
diffusion operator sampled via a heat bath scheme, [isj 
The positive and negative terms are defined by 



B. Methods 

In setting up our Hamiltonian, we use the Born- 
Oppeheimer approximation, a Hartree-Fock norm con- 
serving soft pseudopotential for the He core of carbon 
[60| . and the bare Coulomb potential for hydrogen and 
electron-electron interactions. Our procedure is to start 
with a trial wave function which includes variational pa- 
rameters (see Subsec. lll Al for the forms employed in this 



^r',r = 1/2(Vr'.r±|^R',r|) (10) 

where 

^R'.R - (R'iKon-locallR) , (11) 

and R, R' are all-electron configurations on a quadra- 
ture mesh with one electron rotated around a pseudo 
ion. 44] The breakup corresponds to an effective Hamil- 
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tonian H'^^, defined as 



Tjeff 



off 

R',R 



K + V^iR) 
(R'lK, 



(12) 



-local 



|R) 



if Vr. R < 0, 



witli the modified local potential ^'^^(R) = Vioc(R-) + 
X]r' r that includes the sign flip terms. The 
FN ground state energy of the Hamiltonian in Eq. [T^] 
is a variational upper bound of the original non-local 
Hamiltonian. [45!| Furthermore, DMC stability is im- 
proved substantially due to a softening of the most at- 
tractive parts of the localized pseudopotential. In the 
LA, the highly attractive regions of the localization can 
result in a walker population "blow up" so that the cal- 
culation ends suddenly. Moving the negative part of the 
localization into a diffusion-like term causes the walkers 
to be driven away from such regions. 

The main difference between the DMC Hamiltonian re- 
ported in Eq . [T2l and the LRDMC is the kinetic operator 
K, which is replaced by a discretized K'^ in the LRDMC 
approach, and treated on the same footing as Vnon-iocai- 
K"^ is a linear combination of two discrete operators with 
incommensurate lattice spaces a and a' (a' = i>a, with v 
an irrational number > 1), namely 



(13) 



where A°^p is the discretized Laplacian with mesh a 
and weighting function p (see Refs. [iH and [20l ). and 
r] = 1 + ^a? is a prefactor with the parameter n tun- 
able to improve the efficiency of the diffusion process. 
Working with two incommensurate meshes helps to sam- 
ple densely the continuous space by performing discrete 
moves of length a and a'. The finest hop samples more 
likely regions near atomic centers while the coarser one 
samples more often valence regions, the result being an 
efficient sampling of the overall configuration space. The 
difference between the continuous and discretized local 
kinetic energies is added to y*(R), resulting in a mesh 
dependent potential 



y°(R) = T/'="(R 



(R). 



(14) 



The consequence is a faster convergence of the ener- 
gies in the a — > extrapolation. In spite of the dis- 
cretization of K (Eq. fO)) and the redefinition of V^^ 
(Eq. [M]) , the LRDMC method is equivalent to the con- 
tinuous space FN DMC with Hamiltonian in Eq. [T^l In- 
deed, in the limit of small mesh sizes a and a', the dis- 
cretized Hamiltonian i?" approaches the continuous H . 
The usual DMC Trotter breakup results in a time step 
error while the LRDMC paradigm results in a space step 
error, but both share the same upper bound property in 
the zero-time-step zero-lattice-space limit and converge 
to the same projected FN energy, [isj 

Our SJ calculations were done using continuous space 
DMC with QMCPACK.^ This code provides many fea- 
tures that make it easy to work with SJ wave functions. 



The LRDMC method, available in the TurboRVB,[47| 
has been applied to the JAGP wave function after a 
full optimization of its parameters. We used two op- 
timization procedures. For the SJ work we employed 
the method of conjugate gradients (CG) introduced by 
Hestenes and Stiefel[40| in 1952. This is a first-derivative 
method that finds the minimum of a cost function (in our 
case a linear combination of the variance and the energy) , 
in a number of steps significantly smaller than the stan- 
dard steepest descent method, because for a quadratic 
cost function it converges in a finite number of iterations, 
at most equal to the dimension of the vector space. [48ll49j 
We optimized 10 parameters of the Jastrow functions but 
used the same VTZ basis set at all separations. While 
systems are generally not quadratic, they are so near sta- 
tionary points making the CG approach widely applica- 
ble. Besides this, the method uses little memory, each 
iteration is equal time, and, while not guaranteed, the 
global minimum is usually found. However, the statisti- 
cal noise inherent in the QMC framework limits the appli- 
cability of our CG implementation to systems involving 
not too many parameters, such as our SJ optimization. 

The JAGP optimization, on the other hand, involves 
a large number (^^ 1000) of parameters, mainly coming 
from the (Eq. ^ and g^^ (Eq. O matrices in the 
AGP and Jastrow geminal expansions over the atomic 
basis set. Therefore, an optimization technique robust 
under stochastic conditions is required. The stochas- 
tic reconfiguration (SR) method recently introduced by 
one of us (S.SO 14 1| in conjunction with subsequent 
improvements [l^ l20l. l2ll. [2^ has been shown to be very 
efficient to minimize the variational energy. This is a 
first-derivative algorithm that moves the parameters ak 
at each iteration according to the update 



where 



5ak 



1 

/ . ''k,k' 
fc' 



(15) 



(16) 



Here, k indexes the parameters, fk = —dE/dak are the 
generalized forces, Sk,k' are the SR matrix elements, and 
7 scales the length of the move. The SR matrix is only 
required to be positive definite in order to lower the en- 
ergy. However, to make the convergence more efficient, 
we define the SR matrix as0, IH [2i, |4l| 

sfc,fe' = (1 + 6k,k' e) {{OkOk') - {Ok) {Ok')) (17) 

where Ok ~ da^. In | (RI^'q) | and e is a cutoff chosen 
large enough to guarantee the SR matrix remains well 
conditioned. While the SR matrix is guaranteed to be 
positive definite even for a finite Monte Carlo sample, 
too small eigenvalues can result in an amplification of 
errors coming from the forces fk through Eq. [TBI Set- 
ting e to a finite but small value (usually ~ 10"'^ smaller 
than the largest eigenvalue) makes the optimization much 
more stable. 20] The 7 factor in Eq. (THlcan be tuned to 
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speed up the convergence of standard SR by allowing 
for a better estimation of the magnitude of the parame- 
ter changes 6ak- This is done by evaluating the Hessian 
along the s~^f direction and finding the minimum by as- 
suming quadratic curvature. We extend this idea further 
and minimize the energy by evaluating the Hessian at 
each step over the latest n SR directions, with n chosen 
to maximize the efficiency, as explained in Ref. [20I . 



III. RESULTS 



is supposed to be inert, while in the valence region the 
dependence on the trial wave function is given just by its 
nodes. 

The good agreement between the VMC and DMC 
JAGP results, presented in Subsec flll A) highlights that 
the basis set superposition bias is not relevant (smaller 
than the statistical error of ~ 0.2 mHa) for the fully op- 
timized basis set used in the JAGP wave function, while 
the agreement between the projected SJ and JAGP en- 
ergies, shown in Subsec. IHIBi suggests that the FN bias 
is negligible. 



In this section we present results for hydrogen-benzene 
binding where the hydrogen molecule is oriented along 
the Cg symmetry axis of the benzene molecule. Pre- 
vious studies [H, [l3| found this configuration the most 
stable. Here, we do not take into account other possible 
orientations, because our goal is to check the accuracy of 
different QMC wave functions and provide benchmarks 
for the lowest energy configuration. In order to resolve 
its potential energy surface, we consider the system at 
different molecular center-of-mass separations R. In our 
QMC calculations we have kept the geometry of each 
molecule fixed and close to its experimental structure. 61 
We checked the effect of relaxing the geometries at the 
MP2 level and found an energy lowering on the order of 
/iHa, whose effect is completely negligible in this case. 

We emphasize that all our QMC results do not involve 
any corrections and are direct energy differences with the 
largest computed distance {R = 15) taken as reference for 
the zero of energy. We also present results using PBE- 
DFT where we have quantified the BSSE. In those calcu- 
lations we use the VTZ basis with added diffuse functions 
as described in Sec. Ill Al The BSSE corrected Morse fit 
binding energy and bond length are 0.79 mHa and 6.45 
Bohr. The uncorrected binding energy and bond length 
are 1.18 mHa and 6.25 Bohr. Here, the BSSE is 0.39 
mHa, roughly half the binding. In the QMC framework, 
we do not expect the basis incompleteness error to be 
so important. Although a basis error is unavoidable at 
the VMC level since we use a finite basis set, it is alle- 
viated by the fact that we fully optimize the AGP and 
Jastrow bases along with all exponents at each R. On the 
other hand, the DMC method for local potentials has no 
basis error, because the stochastic projection to the FN 
ground state uses the position representation, which is a 
complete basis set in configurational space, and the only 
bias comes from the FN error. However, in the case of 
non-local pseudopotentials, the localization makes the ef- 
fective DMC Hamiltonian to depend also on the shape of 
the trial wave function (locality error) , and so a residual 
error due to incomplete basis could come through. How- 
ever, the DMC algorithms used in this work are usually 
less sensitive to the locality error, since they are based 
on a partial localization of the non-local potential, as de- 
scribed in Sec. IIIBI Moreover, the locality error, and its 
associated finite basis error, should go away in the energy 
differences, because it affects only the core region, which 



A. Jastrow correlated Antisymmetric Geminal 
Power 



We optimized the variational wave function described 
in Sec. Ill Al bv means of the most recent version of the SR 
energy minimization with Hessian acceleration, f20|] de- 
scribed in Sec. Ill Bl and implemented in TurboRVB.[47j 
The Hamiltonian includes soft pseudopotentials for 
carbon. [131 Although the basis set used here is quite com- 
pact, it turns out that the variational energies are very 
accurate, as we optimize also the exponents of both the 
determinantal and Jastrow part. For instance, the basis 
set for the hydrogen molecule is a {2s2p)/[lslp] Gaus- 
sian in the AGP expansion, while it is an uncontracted 
(Islp) Gaussian plus a constant in the Jastrow geminal 
(a constant generates additional one body terms when 
multiplied by other orbitals Xbm in Eq- [9]) . In spite of 
this small basis set, the variational energy of an isolated 
H2 molecule is —1.174077(29), very close the exact re- 
sult (-1. 174475). [51| The second Gaussian in the s andp 
contractions of the hydrogen AGP is fairly diffuse, their 
exponents ranging from 0.05 to 0.1, as the distance R 
between the benzene molecule and the hydrogen dimer 
shrinks from 15 to 6 Bohr. 

The basis set of the benzene is slightly larger for the 
carbon sites {{6s6p) /[2s2p] in the AGP part, uncon- 
tracted (3s2p) in the three-body Jastrow), while for its 
hydrogen constituents we used just a single s Gaussian 
both in the AGP and Jastrow geminals, since they are not 
supposed to play a key role in the interaction between the 
hydrogen molecule and the benzene ring. This fully opti- 
mized basis set included in the JAGP wave function gives 
a quite good variational energy for aromatic rings. [23| 

We found that the inclusion of the diffuse orbitals 
in the basis set of the hydrogen molecule is crucial for 
the hydrogen-benzene binding, both at the VMC and 
LRDMC level. On the other hand, some Gaussians re- 
lated to the contracted p orbital of the benzene ring be- 
come more delocalized in the binding region. This is rea- 
sonable, because the interaction is supposedly driven by 
the resonance between the carbon pz and molecular hy- 
drogen s components of the total wave function. There- 
fore, the minimal basis set should include diffuse orbitals 
on both sides. We would like to stress that the exten- 
sion of those diffuse orbitals is not determined a priori, 
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but is found by optimizing the wave function with the 
necessary variational freedom. 

After a full optimization of the variational wave func- 
tion at several distances (i? — 5, 5.5, 6, 7, 8, 10, 15) we 
carried out VMC and LRDMC simulations to study 
the properties of the system, in terms of energetics and 
charge density distribution. The LRDMC kinetic param- 
eter in Eq. [13] which optimizes the lattice space extrap- 
olation is = 3.2, that allows one to work with a quite 
large (and highly efficient) mesh size (a = 0.25 a.u.). 
Properly setting the parameters of the LRDMC effective 
Hamiltonian is crucial in order to speed up the simu- 
lation, and so be able to resolve the small binding en- 
ergy of this system. To check the convergence of our 
LRDMC energies with respect to the mesh size, we com- 
puted = 6)-E{R = 15) fora = 0.125, 0.25, and 0.5, 
as reported in Tab. [J It is apparent that the energy dif- 
ferences are converged within the error bar of 0.25 mHa 
in the lattice space range taken into account. It is there- 
fore accurate to work with a ~ 0.25. 



TABLE I: LRDMC binding energy {E{R = 6) - E{R = 15)) 
dependence on mesh size a. The energies are reported in mHa, 
the lengths are in Bohr. 



a 


-S^biiiding 


0.125 


1.53(24) 


0.25 


1.57(19) 


0.5 


2.07(23) 



6.33(15) Bohr for the equilibrium distance, and 1.53(12) 
mHa for the binding energy, as reported in Tab. [Ill 



TABLE IL Fitting parameters of the Morse function (see 
Eq. [IHl) which minimize the of the JAGP-LRDMC and 
SJ-DMC data sets. Their error is computed by means of a 
Bayesian analysis [E^j based on the statistical distribution of 
the FN energy points. The energies are reported in mHa, the 
lengths are in Bohr. 





JAGP 


SJ 


a 


0.56(7) 


0.66(9) 


£b 


1.53(12) 


1.43(16) 


Rq 


6.33(15) 


6.31(21) 
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The results of our calculations of the VMC and 
LRDMC dispersion curves are presented in Fig.[T^, which 
shows the energy as a function of distance R relative to 
the value at R = 15 for each of the methods. There is 
excellent agreement between the two curves, with a dif- 
ference that is less than 0.18 mHa for most points. Of 
course, the diffusion calculation leads to a lower total en- 
ergy than the variational calculation in every case, but 
the agreement of the two methods for the energy differ- 
ence supports the idea that our results are accurate and 
the calculated binding energy is close to the exact value. 

In order to extract the values for the equilibrium 
distance i?o and the binding energy E^, we fitted our 
LRDMC points with the Morse function: 



(18) 



FIG. 1: (Color online) Dispersion energy of the hydrogen- 
benzene bond calculated using various QMC methods for the 
case where the hydrogen molecule is oriented vertically with 
respect to the benzene plane. R is the distance between 
the center of mass of the two fragments. The reference for 
the zero energy difference is taken at i? = 15. The upper 
plot compares the variational and the diffusion results using 
the correlated geminal wave function, labeled JAGP-VMC 
and JAGP-LRDMC. The lower part of the figure compares 
the best diffusion results using two types of trial functions, 
the JAGP (the same as in the upper figure) and the Slater- 
Jastrow function labeled as SJ-DMC. The best Morse fits of 
the diffusion data for the two wave functions are also plotted 
as continuous curves. The close agreement of all three results 
is strong evidence that the binding curve is accurate and the 
analytic JAGP variation function (defined in Eqs. [5H9} is a 
reliable representation of the fully correlated many-body va- 
lence wave function. 



where a is related to the zero point motion of the effec- 
tive one dimensional potential V{R), and E^o is chosen 
to be E{R — 15), i.e. the zero of energy. This choice is 
motivated by the fact that the overlap of the wave func- 
tion in between the two fragments is negligible at that 
distance. Beyond that point the variation of V{R) up to 
infinity is much smaller than the statistical accuracy of 
our points. We estimated the error on the fitting param- 
eters by carrying out a Bayesian analysis of the fit, in a 
way similar to what described in Ref. [s^ . Our result is 



B. Slater- Jastrow Trial Function 

At this point, it is interesting to make a comparison 
with a simple SJ wave function to determine whether 
the use of the JAGP is necessary to get the correct dis- 
persion energy out of the FN projection. We generated 
the Slater part by performing DFT calculations with the 
PBE functional using the VTZ basis ^] with diffuse func- 
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tions from the aug-cc-pVTZ basis. [lOl AH the DFT cal- 
culations were done using Gaussian03.[53l| We chose to 
use a very simple Jastrow factor because our goal was to 
improve DMC efficiency as opposed to obtaining a well 
converged binding curve at the VMC level. Therefore 
our Jastrow includes only one- and two-body contribu- 
tions with three non-cusp terms in the expansion (see 
Eqs. [3]and[3]). Furthermore, a single cusp term is added 
in each of the hydrogen and electron-electron Jastrow 
factors. No cusp is included for carbon due to the use of 
a soft Dseudopotential. j50| 
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-1.0 





LRDMC Morse fit 
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DFT-PBE 
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TABLE III: DMC binding energy {E{R = 6)-E{R = 15)) de- 
pendence on time step r. The energy extrapolated for r ^ 
is within one error bar from the point at r = 0.01. Therefore, 
we chose r = 0.01 as the time step for all our DMC simula- 
tions. The energies are reported in mHa, the time steps are 
in Ha-^ 



r 


Eb 


0.01 


1.38(19) 


0.02 


0.93(19) 


0.04 


0.64(15) 



All the Monte Carlo calculations with the SJ wave 
function were done using QMCPACK.fZ^ The Jastrow 
factor was optimized within the VMC framework us- 
ing the conjugate gradient method. pl| as explained in 
Sec. IIIBI While the SJ variational energy is quite poor, 
its quality is not directly reflective of the DMC energy, 
which depends only on the nodes of the trial wave func- 
tion. In this case, indeed, we found that the DFT nodes 
are very good by carrying out DMC simulations with the 
non-local scheme described in Sec. Ill Bl Our projection 
was done in time steps of r = 0.01 which we found to 
be converged as reported in Tab. IIIII Remarkably, the 
DMC-SJ energies are in very good agreement with the 
LRDMC- JAGP data points (see Fig.flh). Indeed, the SJ 
fitting parameters of the Morse dispersion curve fEg. fTS)) . 
such as binding energy, equilibrium distance, and curva- 
ture, differ from the JAGP ones by less than one error 
bar (Tab. [n|. This consistency between different trial 
wave functions signals that the FN bias is negligible and 
the results are well converged. Moreover, in addition to 
the nodes of the PEE wave function being good, the PBE 
binding energy is underestimated only by a factor 2 with 
respect to our best value. It is notable that the PBE 
functional performs quite well, even though it does not 
include any VdW contribution. In the case of a pure 
VdW bond, the PBE result should be much poorer, as 
already pointed out by Hamel and Cote. [31 This is sug- 
gestive of a more complex binding mechanism which goes 
beyond the standard physisorption. We will focus on this 
point in Sec. [V] 



FIG. 2: (Color online) Results of different theoretical descrip- 
tions of hydrogen-benzene binding as a function of intermolec- 
ular distance R where hydrogen is situated perpendicular to 
benzene. The solid black line and associated error bars show 
the JAGP-LRDMC data and Morse fit (zero binding energy 
is taken at i? = 15 Bohr). The green curve shows the PBE- 
DFT counterpoise corrected result using the VTZ basis [29^ 
supplemented by the diffuse functions from the aug-cc-pVTZ 
basis. ^] The dotted blue curve (shallowest) is the empirical 
potential devised by Crowell and Brown !l5| that takes into 
account the bond asymmetry of the sp^ hybridized carbon 
atom. The dotted red line is the empirical potential by Mat- 
tera et aZ. [3| that seeks to reproduce the hydrogen bound 
states over graphite by a much simpler model. The Mattera 
potential does not take into account bond asymmetry. 



IV. COMPARISON TO OTHER WORK 

The hydrogen-benzene system has been the subject of 
several theoretical works, whereas to our knowledge no 
direct study of this system has been carried out on the 
experimental side. Hydrogen adsorbed on metal-organic 
frameworks (MOF), where benzene- like structures serve 
as ligands, has been studied by Rosi et atHBl who per- 
formed inelastic neutron scattering (INS) measurements. 
The INS data could be related to the rotational states 
of hydrogen adsorbed over benzene. However, the bind- 
ing sites in the MOF structure are not known with cer- 
tainty, and thus it is hard to find a one-to-one correspon- 
dence between the experiment and the isolated hydrogen- 
benzene compound. 

Given the lack of direct experimental data for this sys- 
tem, we compare our results with those from empirical 
models that are often used to estimate complex system 
properties, such as the hydrogen storage capabilities of 
carbon nanotubes and fuUerene nanocages.|54i |55| Here 
we consider two empirical models, both derived from ex- 
periments of hydrogen molecules scattered on graphite 
surface, carried out by Mattera et aZ..[l3| To reproduce 
their data, they proposed a simple model interaction be- 
tween the carbon atoms and the hydrogen dimer which 
depends only on the distance from the graphite layers 
by assuming lateral average. This model was improved 
later by Crowell and Brown, [is'l who constructed an em- 
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pirical potential based not only on the experimental scat- 
tering data but also on the polarization constants built 
in the VdW (6,12) potential. Their model assumes both 
a radial and angular dependence, which takes into ac- 
count the sp^ hybridization asymmetry of carbon atoms 
in graphitic and aromatic compounds. We applied these 
potentials to the hydrogen-benzene system by summing 
the terms for the 6 carbons taking into account distance 
and, for the Crowell potential, the angle the hydrogen- 
carbon interaction makes with the benzene Cg axis. Both 
empirical potentials significantly underbind the system, 
roughly by factors of 3 and 2 respectively when com- 
pared to the J AGP LRDMC results (see Fig. [2]). More 
precisely, Mattera's interaction gives a binding energy of 
0.86 mHa at 5.6 Bohr, while Crowell's gives a minimum 
of 0.54 mHa at 6.2 Bohr. 

Hamel and C6te[l^ calculated the dispersion curves 
using DFT with the local density and generalized gradi- 
ent approximations (LDA and GGA) where the GGA is 
implemented in the PBE density functional. [27l l28j Their 
calculations used a plane wave basis with a 60 Ha cutoff. 
They found that the DFT-LDA gives the strongest bind- 
ing (3.30 mHa), while the DFT-PBE binding is much 
weaker (0.69 mHa). This is consistent with the gen- 
eral overbinding of LDA and underbinding of PBE. It 
is also well known that DFT is not a favorable method 
for systems where van der Waals forces play an impor- 
tant role; [ill in those cases, MP2 and CCSD(T) can 
be applied with more reliability. Hamel and Cote also 
calculated binding curves using those theories. They 
found MP2/6-311-|-G(2df,2p) binding of 1.58 mHa and 
CCSD(T)/6-31+G(d,p) binding of 0.65 mHa. 

Perhaps the most careful and accurate MP2 and 
CCSD(T) calculations were done by Hiibner et aZ.jl3] In 
order to resolve the weak interaction between hydrogen 
and benzene, high accuracy is required, and so a large 
basis set is needed to reduce both basis set superposition 
and incompleteness errors which are a significant frac- 
tions of the binding energy (the BSSE was found to be 
as much as ^ 25% of the final estimated binding). On 
the other hand, the use of a larger basis set is limited by 
a poorer scaling of the calculations, particularly at the 
CCSD(T) level of theory, which is the most expensive. In 
their work, Hiibner et al. optimized the binding distance 
using MP2 with the TZVPP basis. They found a center- 
of-mass distance of 5.80 Bohr and a binding energy of 
1.47 mHa. This geometry was then used for further MP2 
and CCSD(T) calculations. The CCSD(T) method with 
the same TZVPP basis gives 1.17 mHa, while the MP2 
theory was pushed up to a aug-cc-pVQZ' basis to give 
a binding of 1.83 mHa, a significant increase from the 
TZVPP basis. At this point, it is possible to estimate 
the true binding energy by correcting the best MP2 en- 
ergy with the CCSD(T)-MP2 difference obtained at the 
TZVPP level. This gives a value of 1.5 mHa, remark- 
ably close to the JAGP LRDMC binding of 1.53 ± 0.12 
mHa, found in this work. 



V. ANALYSIS OF THE BONDING 

In order to investigate more deeply the physics of hy- 
drogen adsorbed on benzene, we study the induced dif- 
ference in electronic density at the equilibrium bond dis- 
tance with respect to the separated fragments. For this 
study we compare our best DMC results to the density 
functional calculation using the PBE functional. The 
QMC densities are calculated from the optimized cor- 
related geminal (JAGP) as a mixed estimator, which is 
an accurate representation of the DMC results since the 
diffusion calculation leads to only small changes (within 
the error bar) from the VMC density. The contour plot 
in Fig. [3] shows the difference in the calculated electron 
density at the separation i? = 6 Bohr. Here, the elec- 
tron density of the isolated molecules has been subtracted 
from the combined system so that the change in charge 
distribution due to bonding is apparent. In this figure 
the benzene ring lies in the xy plane at z = and the 
hydrogen molecule is oriented along the z axis, with its 
center of mass at z = 6 Bohr. The two dimensional plot 
in the yz plane is generated by integrating the density 
distribution over the x coordinate. As one can see, the 
hydrogen molecule is polarized by the electronic repul- 
sion with the benzene cloud, which pushes the electrons 
to the opposite side of the molecule, leading to a static 
dipole moment on the hydrogen. On the other hand, the 
density redistribution in the benzene is non trivial, and 
shows patches of charge accumulation and depletion. To 
catch the net effect of this redistribution, we integrated 
the density also over the y coordinate, and obtained an 
effective linear density profile, plotted in Fig. 21 Here, it 
is apparent that the overall effect on the benzene is the 
formation of another effective dipole moment, oriented 
to the same direction as the static dipole moment on the 
hydrogen molecule, which lowers the electrostatic energy. 
Notice that in Fig.|3]we have plotted separately the VMC 
and the LRDMC mixed estimate of the densities. The 
close agreement supports our conclusion the VMC wave 
function is very accurate not only for the energy but also 
for other properties such as the density. 

At large distances the attractive interaction is due to 
VdW dispersive forces, which is included in the Monte 
Carlo calculations. At short distances the interaction is 
repulsive due to overlap of the closed shells, which would 
lead to density displaced outward on both the hydrogen 
and benzene, i.e. opposite dipoles on the two molecules. 
However, Figs. [3] and H] show that the hydrogen-benzene 
bond is not a pure VdW interaction, since in the binding 
region also electrostatic effects come in with the onset of 
dipolar interactions that lower the charge repulsion. For 
comparison, density differences calculated using the PBE 
density functional are also shown in Figs. [3] and [Hat the 
separation R = 6 Bohr. Of course, the PBE functional 
does not include VdW interactions so that the binding 
decreases too rapidly at large distance as shown in Fig. [21 
Nevertheless, near the equilibrium distance the density is 
similar to the QMC result but with smaller magnitude 
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of the change in density, which is consistent with the 
fact that the PBE functional underbinds the system. It 
is well known that GGA functionals like PBE tend to 
underbind because they favor systems with larger gra- 
dients, whereas LDA tends to overbind molecules and 
solids since it favors more homo gene ous systems. [49| Re- 
cent work by Langreth et al. 5^, 53| has led to improved 
functionals including van der Waals interactions; how- 
ever, they have not been considered here. 
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FIG. 3: (Color online) Contour plots of the difference in 
projected electronic charge per unit area between hydrogen- 
benzene separated by 6 Bohr and the isolated hydrogen 
and benzene using JAGP-LRDMC and PBE-DFT. This plot 
shows how the charge per unit area changes as hydrogen and 
benzene interact, the x-axis have been integrated over so that 
the charge per unit area has been projected into the yz-plane. 
(Left) The areal charge density difference is a mixed estimate 
of LRDMC calculations with a JAGP trial wave function. 
(Right) computation is done within the PBE-DFT framework 
using the VTZ basis [2^ supplemented by diffuse functions 
from the aug-cc-pVTZ basis, [s^ The combined hydrogen- 
benzene system basis was also used in the isolated benzene 
and hydrogen calculations. 



VI. SUMMARY AND CONCLUSIONS 

We have presented VMC and DMC results for the ad- 
sorption of hydrogen on a benzene ring and compared 
them with previous work. We used two types of varia- 
tional correlated wave functions, a S J function with DFT- 
PBE optimized single-body orbitals and a JAGP function 
fully optimized at the VMC level by means of the SR en- 
ergy minimization. In this work, we have shown strong 
evidence that our results are very accurate since we have 
found essentially the same results in three independent 
QMC calculations: one J AGP- VMC variational simula- 
tion with no FN error, and two DMC simulations based 
on different trial wave functions (JAGP and SJ) with pos- 
sibly different nodes and no basis set errors. The agree- 
ment among our three calculations is within ~ 0.2 mHa, 
which is mainly due to statistical accuracy on the QMC 
energies, and gives an upper bound for the magnitude of 
underlying errors, such as the basis set incompleteness 
and the FN bias. 
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FIG. 4: (Color online) Plot of the difference in linear elec- 
tronic charge density between hydrogen-benzene separated by 
6 Bohr and the isolated hydrogen and benzene using three 
theories. This plot shows how the linear charge density dis- 
tribution changes as hydrogen and benzene interact, the x- 
and y-axes have been integrated over so that only the z-axis 
is shown. The positions of the hydrogen and benzene are in- 
dicated in the graph. The solid red data with error bars show 
the induced charge changes using the analytic JAGP wave 
function at the VMC level. The dotted blue data with er- 
ror bars show the mixed estimate of the density, given by the 
LRDMC projection of the JAGP trial wave function. The dot- 
ted green line without error bars shows the PBE-DFT charge 
density difference using the VTZ basis [2^ supplemented by 
the diffuse functions from the aug-cc-pVTZ basis. [30l| The 
combined hydrogen-benzene system basis was also used in the 
isolated benzene and hydrogen calculations. 



Our best estimate for the binding energy is 1.53(12) 
mHa at an equilibrium distance of 6.33(15) Bohr, ob- 
tained by using the LRDMC method with the nodes of 
the JAGP wave function. Our result agrees well with 
the conclusion of Hiibner et a/. [13] who used MP2 and 
CCSD(T) methods, and estimated the binding to be 
~ 1.5 mHa based on extrapolation which accounts for 
basis set and level of theory. The resulting binding en- 
ergy is 2-3 times larger than those given by empirical 
potentials [13, [3 ^^'^ DFT-PBE calculations which are 
often employed in more complex systems, suggesting that 
their results could be substantially affected by this lack of 
accuracy. It would be interesting to extend the present 
work, by studying the transferability of such empirical 
potentials on other aromatic and graphitic structures. 

We proved that the JAGP wave function provides a 
very accurate dispersion curve for this system already 
at the variational level. This result is remarkable, be- 
cause we were able to derive a compact analytic form 
which can be used for accurate determination of proper- 
ties other than the energy by means of the VMC method 
with no sign problem. The JAGP wave function cap- 
tures the resonating valence bonds of benzene in its gem- 
inal construction as well as the van der Waals interaction 
through many-body correlations in the Jastrow factor, as 
shown in previous work on benzene dimer.*2Cf| The basis 
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for both the AGP and Jastrow geminals is of compact 
Gaussian form that does not go beyond p-orbitals, but 
includes diffuse orbitals with optimized exponents. 

By means of the DMC method, we also studied the 
hydrogen-benzene problem using a more conventional SJ 
wave function. The single-body orbitals included in the 
Slater determinant were derived using the GGA PBE 
density functional with a VTZ Gaussian basis [l^ mod- 
ified to include diffuse functions from the aug-cc-pVTZ 
basis, [soj as discussed in Subsec. Hi Al which is essential 
for an unbiased DMC sampling of long-range VdW ef- 
fects. The Slater basis set is roughly 4 times larger than 
its JAGP counterpart whereas the Jastrow factor is of 
minimal form, satisfying cusp conditions and improving 
computational scaling. Our findings suggest that for this 
particular problem the geminal form is not essential to 
get an accurate DMC energy, and can be replaced with 
a Slater determinant and DFT optimized orbitals in the 
DMC calculations. While the JAGP uses a more com- 
pact basis, the SR optimization involves a large number 
of parameters coming from the Jastrow and AGP gem- 
inals expanded on atomic orbitals. This makes the SJ 
wave function with DFT-PBE single-body orbitals more 
desirable for DMC calculations in larger related systems. 

Finally, we examined how the electronic density of 
the isolated molecules changes in the bond region. The 
change in density, displayed in Figs. [3] and [H shows that 
near the equilibrium distance there is the formation of 
static dipoles that can lower the electrostatic energy, in- 
dicating a bonding mechanism beyond VdW. Density 
functional calculations using the PBE functional lead to 
similar density profiles but with smaller magnitude, in 
agreement with the well known underbinding tendency 
of that functional. This means that the interaction be- 
tween hydrogen and benzene is not a pure VdW effect. 



since it can be partially captured by a DFT-PBE formal- 
ism which does not include dispersive interactions. This 
also clarifies why the DFT-PBE nodes of the SJ wave 
functions are very good, and equivalent to the JAGP 
nodes to predict the correct binding energy at the DMC 
level. 

To conclude, we have reported on a detailed analy- 
sis of the hydrogen adsorption over molecular benzene 
by QMC methods, which are shown to be very accurate 
and reliable to predict the energetics and other physical 
properties of the system. This framework is therefore 
promising to study hydrogen interacting with graphitic 
or other aromatic compounds, particularly important for 
the hydrogen storage problem. 
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